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Abstract 

A new non-perturbative method of solution of the nonlinear Heisenberg 
equations in finite-dimensional subspace is illustrated. The method, being a 
counterpart of the traditional Schrodinger picture method, is based on a finite 
operator expansion into the elementary processes. It provides us with the 
insight into the nonlinear quantal interaction from the different point of view. 
Thus one can investigate the nonlinear system in both pictures of quantum 
mechanics. 

1 Introduction 

The use of laws of quantum mechanics in the description of nonlinear systems con- 
fronts us with the qualitatively new difficulties. Namely, to investigate their dynam- 
ics in the Heisenberg picture we have to solve the nonlinear operator equations, a 
task which is highly nontrivial even for the simplest systems. The difficulties are also 
encountered with in the Schrodinger picture once we try to solve the Schrodinger 
equation explicitly fl. Since some nonlinear systems solvable analytically in clas- 
sical domain become insoluble when are quantized one can suppose that they are 
simultaneous influence of intrinsic stochastic effects, originating from the incom- 
patibility of some observables, and nonlinearity which make the behaviour of such 
systems very complex and thus difficult to describe analytically. 

The time evolution of quantum systems can be studied with the help of widely 
used Schrodinger picture method based on the integration of set of linear differential 
equations for components of a state vector in the Fock basis Q , ||| . Unfortunately, 
the expansion into the Fock-state basis can be infinite for some states, e.g. for 
coherent state, yielding the infinite set of these equations. Because it is practically 
impossible to solve the infinite system of the equations, the method provides us 
with exact solutions only for states from some finite-dimensional subspace of the 
Hilbert state space. On the other hand, it is advantageous sometimes to calculate 
the evolution of particular observables in the framework of Heisenberg picture. The 
motivation of the present paper is to find the operator analogue of the Schrodinger 
picture method in the Heisenberg picture and to show the equivalence and deep 
relationship between them. 

As an illustrative example we consider here the simple system composed of 
two harmonic oscillators which oscillate with frequencies u> and 2u> and which are 
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described by the annihilation (creation) operators ai(a|) and 02 (aj) obeying the 
standard boson type commutation rules 

[oi,at] = Sij, [ai,Oj] = [Si.a]] = 0, i,j = 1,2. (1) 

Let the two oscillators interact nonlinearly according to the following interaction 
Hamiltonian (the second-harmonic generation process Q ) 

H = -hTa\ 2 a 2 + h.c. , (2) 

where T denotes the nonlinear coupling constant; the symbol fi is reduced Planck 
constant and h.c. stands for the Hermitian conjugate term. Here and in the follow- 
ing we assume that the free evolution was eliminated by the appropriate unitary 
transformation. 

Employing the commutation rules (Q) one can directly prove the existence of the 
following integral of motion 

7V = ni+2n 2 , (3) 

corresponding to the total energy of the system. Here the photon number operator 
hj of the jth oscillator , j = 1, 2 has been introduced. The eigenvectors of the 
integral of motion (|^) then provide us with the natural orthonormal and complete 
basis in which the expressions have simple form. They are easy to find and have 
the form 

riVi 

{\N-2l,l),l = 0,l,..., [— J,iV = 0,l,...} (4) 
with orthonormality condition 

(N-2l,l\M-2k,k)=6 NM 5 lk (5) 
and the resolution of unity operator 



J2J2\N-2l,l)(N-2l,l\ = l, (6) 

jv=o 1=0 

where Jia) is the Fock state having energy htoni + 2hton2', N is the eigenvalue 
of (H), [N/2] represents the greatest integer less or equal to N/2 and Sik is the 
Kronecker symbol. The Hilbert state space of our system can then be expressed as 
a direct sum 

00 

h = ® n(N) ( ? ) 

N=0 

of the invariant [iV/2] + 1-dimensional subspaces HS N ' spanned on the basis vec- 
tors \N — 21, 1), I = 0, 1, [N/2] corresponding to the fixed eigenvalue N. Using 
the standard properties of the annihilation and creation operators of the harmonic 
oscillator 

ai\rii) = y/nl\ni-l), a\\rii) = y/m + l\m + 1), i=l,2 (8) 

and employing the condition (^|) one can show that the Hamiltonian (^) is repre- 
sented by the following block diagonal matrix 



(N -2l,l\H\M -2k,k) = -h\Ty/(l + 1)[N - 21) (N -21- l)5 k>l+1 

+Vy/l(N -21 + 2)(N -21 + 5 NM , (9) 

where the symbol '*' represents the complex conjugation. 
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2 Schrodinger picture 



Let us first recall the results obtained with the help of the Schrodiger picture method 
when applied to our system. As is well-known the time evolution of the state vector 
is governed by the Schrodinger equation 

mj t \m) = mm), (io) 

where Hamiltonian H is given in (|^). Rewriting ( |Tc| ) into the basis (|4|) with the 
help of (^) and (||) we successively arrive at the infinite number of sets of linear 
differential equations 

- iT^/(l + l)(N-2l)(N-2l-l)C N , l+1 



+iT*y/l(N -21 + 2){N -21 + ljC^-i (11) 

for components C N j = {N — 2l,l\ip(t)), where N = 0,1,... and I = 0, 1, [N/2]. 
Assuming, however, the initial state to be from the finite-dimensional subspace 

K 

H K = ® niN) > (I 2 ) 

N=0 

it is sufficient to solve only K + 1 such sets labelled by eigenvalues N = 0, 1, K 
each of them with [AT/21 + 1 equations. Particularly, for states belonging to the 
subspaces 7i 2 , the set ( |ll|) is of the form 

and can be solved analytically. An interesting result is obtained assuming the system 
to be in the state |0, 1) at the beginning of the interaction. The initial conditions for 
the set (||) are then C ,o(°) = Ci.o(0) = C 2 . (0) = 0, C 2 .i(0) = 1 and the solution 
of @ reads 

2 [-f 1 

W*)> = EE C N +{t)\N - 21, l) = i^- sin(V2|r|i)|2, 0) + cos(V2|r|t)|0, 1). (14) 

N=0 1=0 ' ' 

Hence we obtain the following expressions for the mean number of energy quanta 
in oscillators 1 and 2 in state ( |l4[ ) 

(ni(t)) - 2sin 2 (V2|r|<), (h 2 {t)) = cos 2 (\/2|r|i). (15) 

This non-classical oscillatory behaviour can be interpreted from the point of view 
of Schrodiger picture as being a manifestation of quantum interference effect (|l4j) . 

3 Heisenberg picture 

In this picture the operators d\ and d 2 for the system of interest evolve according 
to the Heisenberg equations of motion 

ih^ = [a j ,H] ) j = 1,2, (16) 
which after substitution (0) into (|l^) and application (|l]) read 

— =2iTa[a 2l — = iT a x . 17 
dt dt 
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It is also well-known that the operators a±(t) and a 2 (i) can be equivalently expressed 
as follows 



a 3 (t) = exp(^Ht)a 3 (0)exp(-^Ht) = £ j = 1,2, (18) 

fc=0 



where the exponential operators have been expanded and (|16|) has been used re- 
peatedly. 

The power series on the right hand side (R.H.S.) of ( |l8| ) is a perturbative so- 
lution of the equations ([l^) and provides us with two important informations. As 
becomes clear from the following it is advantageous to work in the normal ordering 
of the operators in which all creation operators stand to the left from all annihila- 
tion operators. First, the operator part of the solution (|lj), given by derivatives 
of operators a\(i) and a 2 (t) at t = 0, cannot contain products of operators other 
than those of leading to the annihilation of one energy quantum from the corre- 
sponding oscillator. This can be proved by deriving the Heisenberg equations of 
motion ( |i~7| ) and using consequently the commutators (|l]) to obtain the normally 
ordered expressions. From now any such product of operators at t = is called a 
process in the corresponding oscillator and the number of operators in the product 
is called an order of the process. Secondly, calculating the perturbative solution 



(18) to the sufficiently high order and rearranging its terms appropriately, one can 
see that the solution is of the form of finite sum of the processes multiplied by 
various polynomials in t, which can constitute the first few terms of power series of 
well-known functions (it can be verified at least for the first few processes). This 
different point of view to the standard perturbative solution Q is the core of our 
non-perturbative method developed in the following text. Hence one can surmise, 
that going to the infinity in the iterative procedure, the solution of the Heisenberg 
equations of motion ( |l7| ) is of the form of infinite sum of processes multiplied by 
some time dependent functions 



ai 



(t) = + f 1 {t)a\a 2 + f 2 {t)a\al + f 3 {t)a f 2 a 1 a 2 + (19) 



a 2 (t) = gi(t)a 2 + g 2 {t)al+ g 3 (t)a\a 1 a 2 + g4(t)alaj + ... (20) 

where a,j = a,j(0), j = 1,2. The functions fj and gj are called amplitudes of the 
corresponding processes in the following text. 

Substituting ( |l9| ) and (20) into ( |l7| ) and comparing the coefficients related to the 
same process, the amplitudes fj and gj can be determined as solutions of a system 
of ordinary differential equations. For example, the equations for amplitudes /i and 
gi together with the initial conditions read 

j t h{t) = 2£T gi (t), j t9l {t)=iT*h{t), /i(0)=0, «/i(0) = l (21) 
and have the following solutions 

/i(t)=^sin(V2|r|t), Sl (t)=cos(V2|r|t). (22) 

Employing ([l9]) and (|o|) the operators of the number of the energy quanta in the 
oscillators 1 and 2 are of the form 

n^t) = a\ai + 2sm 2 (V2\T\t)ala 2 + (23) 
n 2 (t) = cos 2 (V2|r|f)a3,a 2 + (24) 
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where relations (lj) and ( p2| ) have been used. It is worth noting, that contrary to 
the R.H.S. of (p4) the second term in ( p3[ ) is the quantum contribution originating 
from the commutator [ai,aj] = i. Considering as in the previous section the input 
state to be |0, 1) state and assuming that the next terms represented by dots in ( p3| ) 
and ( p^ ) do not contribute, one obtains for the mean number of the energy quanta 
in the oscillators the expressions 

(ni(t)) = \h(t)\ 2 = 2sin 2 (V2|r|t), (n 2 (t)) = \ 9l (t)\ 2 = co S 2 (V2\T\t), (25) 

which are identical with the results ( |l5| ) obtained by means of the Schrodingcr 
picture method. Notice, that this derivation illustrates not only the mathematical 
equivalence of both methods but also their difference when one tries to distinguish 
between the classical and quantum contributions. 

Although one could look at the method just described as being a satisfactory 
method, let us recall the reader, that its conclusion (|25| ) rests on two crucial as- 
sumptions which were not justified at all. First, we have assumed implicitly, when 
deriving (|2l]), that the higher order processes do not affect the first order ones 
(cosequently we have obtaind the finite set of differential equations for amplitudes 
fi and gi). Secondly, the mean numbers of energy quanta in state |0, 1) given by 



(25) have been derived under the assumption that only explicitly given terms in 



(23) and fl24|) contribute. To show that this is really the case, we have to formalise 



and precise the Heisenberg picture method. This is done in the following section. 



4 General method 

The previous section provides us with an illustrative example, how one can treat 
the system (^|) within the framework of Heisenberg picture on the intuitive basis. 
In the present section we try to justify the intuitive assumptions discussed above 
and to generalize this treatment to the arbitrary finite-dimensional subspace. This 
can be achieved by the suitable parametrization of the problem under discussion. 
To that aim let us rewrite the expansions (|l^) and (^) into the compact forms 

oo 

Si(t) = E fmimalTi&ly^A, (26) 

i,j,k,l=0 
oo 

fi a (t) = £ 9oprs {t){a\)°{alYa\a^ (27) 

o,p,r,s— 

where 

21 + k - 2j - i = 1 (28) 

and 

2s + r - 2p - o = 2 (29) 

holds. There are two facts which can make the convenient parametrization easier 
to find. First, as in the Schrodinger picture we can employ the existence of the 
integral of motion (||). Secondly, the discussion in the previous section indicates 
that the order of the process is of importance. Therefore we put 

N = k + 21, M = i + 2j, m = k + l, R = i + j + k + l, (30) 

for oscillator 1 and similarly 

N = r + 2s, M = o + 2p, m = r + s, R = o + p + r + s, (31) 
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for oscillator 2. Thus each process is parametrized by the parametres N(M) (rep- 
resenting the amount of the annihilated (created) energy in the process), to (the 
total number of the annihilated energy quanta) and R (the order of the process). 
Substituting ( |3C| ) into ( |26| ) and eliminating N by means of (|2^) we arrive at the 
following expansion 

oo M+l M+m 
ai(t) = E E E fMmB.{t) 

x{a) l f R - M - 2m {a\) M+m - R a\ m - M - 1 a^- m+1 . (32) 
In the same way we obtain 

oo M+2 M+m 

fi 2 (*) = E E E 9MmR(t) 

M=0 m=[ M±i]+i R =[M+l] +m 

x(a t 1 ) 2J? - M - 2m (4) M + m - fl a^-M-2 a M- m +2^ (33) 

Substituting ([32|) and (|33| ) into the Heisenberg equations of motion (0), using the 
following equation g 

mm(m.n) , \ / \ 

a?(alT= g '"!" i = l,2, (34) 

which can be proved easily using the commutation rules ([!]), and comparing the 
expressions corresponding to the same process, we obtain the following infinite set 
of differential equations for amplitudes /mto-Ri M — 0, 1, m = [4f] + 1, ...,M+ 1; 
i? = [^i] + m, M + to and ffMm fl, M = 0, 1, to = [Mil] + 1, M + 2; 
#= [^] +m,...,M + m: 

, oo Mi + 1 A/2+2 Mi+mi M 2 +m 2 

|ww - *r x; E E E E 

M ll M 2 =0 mi=[ ^ ]+lTO2=[ ^ ]+1Hl=[ ^ ]+mifl2=[ M^i ]+TO2 

2Ri - Mi - 2mA /2i? 2 - M 2 - 2m 2 \ /Mi + m x - i?i N 
s\ J \ s\ J \ s 2 

2 + TO2 2 j Sl ! S2 !/* fimiKi (t) 5M2m2i?2 (t) ; (35) 

S 2 / 



, 00 Mi+1 M 2 + l Mi+mi M 2 +m 2 

3Af m fl(i) = * r * E E E E E 



X 



2toi - Mi - 1\ /2i? 2 - M 2 - 2m 2 \ /M x - mi + 1 
si / I si A s 2 



x , M 2 + to 2 2 I Sl | S2 !/ Mimifli (t)f M2 m 2 R 2 (*), (36) 



where 



s 2 



i? 2 = R + Ri - 2m - 2toi + 2m 2 , 
si = M — Mi — M 2 — 2to — 2toi + 2to 2 + 2i?i — 1, 
s 2 = Mi + M 2 - M - to - toi + to 2 + R - R 2 + 1 (37) 
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for equation (pqj and 



i? 2 = R - Ri - 2m + 2m 1 + 2m 2 , 

si = M - Ml - M 2 + 2m - 2m a - 2to 2 - 2i? + 2i? 4 + 2R 2 , 

s 2 = Mi + M 2 - M - m + m x + m 2 + i? - #x - #2 (38) 



for equation (p6|). The initial conditions for equations (35) and ( |3q ) are /oxx(O) = 
3oxx(0) = 1, /Mmfl(O) = 9MmR(0} = in all other cases. 

The structure of R.H.S. of (|3g) reveals that there is a nonzero contribution to 
the R.H.S. only if the following inequalities hold simultaneously 

2i?x-Mx-2mx > s x , 2i? 2 - M 2 - 2m 2 > a x , 

Mx + mx - Ex > «2, M 2 + m 2 - i? 2 > s 2 . (39) 

Combining this with one finally obtains 

Mx < M — 1, M 2 < M — 1. (40) 

Analogously, R.H.S. of ( |3^ ) contains nonzero contribution only if 

2mx-Mx-l > sx, 2i? 2 - M 2 - 2to 2 > sx, 
Mx-mx + 1 > s 2 , M 2 + m 2 - i? 2 > s 2 , (41) 



hold simultaneously. Consequently, the R.H.S. of ( |36| ) cannot contain amplitudes 
other than those for which 

Mx<M, M 2 <M + 1. (42) 



From the above inequalities (4C) and (|4^) follows that for fixed M we have only 
finite set of equations (|35|), (36) for amplitudes /jw'm'H'jM' = 0,1,..., M; m' = 



[^-] + l, M' + l; i?' = [MjtlJ+TO^ M'+m' and g M "m"R",M" = 0, 1, M-l; 
m" = i^^ 1 ] + 1, ...,M" + 2; i?" = [^5^] + m", M" + m". In other words, 
the process in the oscillator 1(2) parametrized by M' > M(M" > M — 1) does not 
affect the processes in the oscillator 1(2) for which M' < M(M" < M — 1), as we 
wanted to prove. 

The discussion of the structure of the set of equations (|35|), ( |36| ) can go even 
further. Since the amplitudes fM'm'R',M' < M and gM"m" R." 1 M" < M ~ 1 can 
be calculated solving the set of equations (|35|), ( |3^ ) corresponding to M — 1, in fact 
only the amplitudes f Mm 'R',m' = [4f ] + 1, M+1;R'= [^tl] + to', A/ + to' 
and gM-im"R",m" = [f-] + l, ...,M+1; i?" = [4f]+m", ...,M+m"-l are mutually 
coupled. The amplitudes fwm'R' ,M' < M and gM"m"R" >M" < M — 1 then play 
the role of known coefficients and source terms and the set of differential equations 
corresponding to M is linear. Moreover, taking M and to fixed, substituting ( p7| ) 
into ([59|) and putting Mi = M 2 = M — 1 one obtains the following equality 

to 2 = to. (43) 

Since the same equality can be proved substituting ( ^8| ) into (^lj) and putting Mi = 
M, M 2 = M + 1, one can conclude that only amplitudes fMmR/,R' = [ o" 1 ] + 
to, M + to and gn-imR", R" = [-y] + to, M + to — 1 are coupled. Hence for 
given M and to one has to solve the set of M + 1 differential equations (^5|) and 
)• 

Before going further let us notice that since the infinite series (|32|) and (B3T) with 



(36) 



amplitudes being the solutions of the equations ( |35| ) and (36) satisfy the Heisenberg 
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equations of motion ( p"7j ) identically, the operators &i(t) and a 2 (t) preserve the 
commutation rules (|l|). 

Someone still could object that our method cannot be used in practice since we 
are not able to calculate all the amplitudes. This difficulty is, however, overcome 
if one realizes the following fact. Calculating the matrix element of the process 
corresponding to parameters M, m and R for oscillator 1, 

(JVi - 2/ 1 ,/ 1 |(al) 2fl - M - 2m (4) M + m - fl a2 m - M - 1 af- m+1 |7V2 - 2l 2 ,h) 
I ^i!Z 2 !(iVi-2/ 1 )!(iV2-2^)! 

W-2m + l)\ 

(44) 



y (Nt - 2h +M+ 2m - 2R)\{N 2 - 2l 2 + M - 2m + 1)! 

fiNi-2h+4m,N2-2h+2R+l8li+R+l,h+2m 

y/(h — M — m + R)\{1 2 - M + m - 1)! ' 

where the formulas (H) and d|) have been used, it is evident that it does not vanish 
only if the following inequalities are satisfied simultaneously 

Ni - 2h > 2R-M- 2m, 2h > 2M + 2m - 2R, 

N2-2I2 > 2m-M-l, 2l 2 >2M-2m + 2. (45) 

Hence 

Ni > M, N 2 >M + 1. (46) 

Repeating the same discussion for the same matrix element of the process in oscil- 
lator 2 characterized by the parameters M, m and R one arrives at 

Ni > M, N 2 > M + 2. (47) 

The inequalities ( f46| ) and ( p7^ can be interpreted as follows. Restricting ourselves to 
the finite-dimensional subspace Hk of the whole Hilbert space (|7|) only processes in 
oscillator 1 (2) for which M' = 0, 1, K— 1 (M" = 0, 1, K- 2) are represented by 
nonzero matrix. In other words, the time evolution of the operators a±(t) and a 2 {t) 
on the subspace Hk is known once the amplitudes ju'm'B' , M' — 0, 1, ...,K — 1; 
m' = [^f] + 1,...,M' + 1; Bf = [M^k] + m ',...,M' + m' and g M "m"R",M" = 
0, 1, K — 2; m" = [^±i] + 1, ...,M" + 2; R" = [M^tl] + m " , ...,M" + m" are 
determined. This requires sequential solution of [^] + 1 sets of M = 0,1, ...,K 
differential equations ([35]) and ([36|). Since the series (|3^) and (^3|) are terminated 
naturally when considering only finite-dimensional subspace TLk, a natural question 
arises whether the last-named amplitudes determine not only the evolution of ai (t) 
and a 2 (t) but also the evolution of any operator on the subspace TLk ■ Now we will 
prove that this is really the case. 

It is well known that any operator at time t on the space 7i can be expressed as 
a sum of the following products (a\) 1 {t){a\y {t)a\(t)a 2 (t). Hence, it is sufficient to 
prove the statement for these products only. The commutation rules (|l|) enable us 
to show that 

^ (a t ) 2fl-M- 2m(a t)M +m -Ji a 2m-M-l a M-m+l| iv _ ^ q 

= (N- l)( a t)2K-W-2 TO(a tjJlf+ m -« a 2m-M-l a M- m +l| Ar _ ^ (4g) 

for oscillator 1 and similarly 

iV(al) 2fl - M - 2m (a 2 ) M+m - fl a 2m - M - 2 a 2 vf - m+2 |iV - 21, 1} 

= (N — 2)(a t 1 ) 2J? - M - 2m (4) M+m - J? a 2m - M ' 2 af ~ m+2 \N - 21, 1) (49) 

for oscillator 2. Consqucntly, 

Na,{t)\N -21,1) = {N-j)a j {t)\N-2l,l), j = 1,2, (50) 
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as one can verify, using ( p2| ) and (|33|). From that it follows that aj(t)H.K C Tix-j, 
j = 1,2. Since the first annihilation (creation) operator to the right (left) in the 
matrix elements 

(iVi-2i 1 ,« 1 |(al) i (i)(a|) i (t)aJ(t)a^(t)|iV 2 -2; 2 ,Z 2 ), JVi, 7V 2 = 0, 1, K (51) 

transforms the basis vector to the right (left) into the subspace embeded into TCr, 
the series ( |32| ) and ([53]) for the following annihilation (creation) operators must 
terminate even further than those for the first annihilation (creation) operator. 
Therefore no other amplitude except for those mentioned above can appear in the 



expression (51). This is what we wanted to prove. 

There is one more point connected with the previous discussion which should 
be clarified here. Namely, one could think about the finite series for a\(t) and a 2 (t) 
on the TLk as an approximate operator solutions of ([l7j). The following special 
example disproves the idea. 

Let us consider the following finite series (describing correctly the time evolution 
on the subspace H 2 ) 

oi(t) = fonit)^ + f 112 (t)a\a 2 + f 123 (t)a\al (52) 
a 2 (t) = gou(t)a 2 +g 02 2(t)af. (53) 



The amplitudes in (52) and ( |53| ) are solutions of the set of equations 
|/onW = 0, 

j t hi2(t) = 2ir/ *n(t)ffoii(t), ^S0u(t) = ir7ou(t)/ii2(t), 

^/i 23 (t) = 2ir/ * u (*)fl(B2(*), ^soa2(t) = ir7ou(t)(/oii(*) + /i23(t)) ) 

(54) 

with the initial conditions /on(0) = ffou(0) = 1, /m(0) = /i 2 a(0) = 3022(0) = 0. 
The use of standard methods then yields 

/oii(*) = 1, 

f 112 (t) = i^— sm(V2\T\t), sou(t) = cos(V2|r|t), 

/ ia3 (t) = cos(v^|r|*) - 1, g 22(t) = i-^-sm(V2\r\t). (55) 

v2|l I 



Now, substituting (|5J]) into (^3|) one arrives at 

[a 2 (t),al(t)\ = l + 2sin 2 (V2|r|t)alai ^ i. (56) 

In the same way we can prove that [a\{t) , a\(t)] ^ i. This illustrates that the 
commutation rules are not preserved for the finite series ( |52| ) and (^3|). One can 
expect this, since the series (^) and (|5^) represent a correct solution only on the 
subspace H 2 . Thus only the complete solution ([32]) and (|3^) involving all the 
processes preserves the commutation rules in the whole Hilbert space. 

5 Conclusion 

On the simple example we illustrate how to solve the nonlinear Heisenberg equa- 
tions of motion on the finite-dimensional subspace using the finite expansion of 
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annihilation operators into the sum of the elementary processes. The idea of the 
method is not restricted to this example and provides us with recipe how to treat 
other nonlinear interactions. The time evolution of any operator on the subspacc 
is then governed by a finite number of the c-number differential equations for am- 
plitudes. Due to the hierarchy of the processes the equations split into several sets 
which can be solved step by step. Thus the problem of solution of the q-number 
Heisenberg equations is transformed into the finding of solution of the linear c- 
number differential equations, which can be handled numerically. It provides a nice 
interpretation and deeper insight into what happens in the course of the nonlinear 
quantal interaction in the language of elementary processes. It also enables us to 
identify the non-classical contributions. This instructive interpretation cannot be 
obtained within the framework of the Schrodinger picture. 
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